Erwinia (oligotyping)

Load packages, paths, functions

# Load main packages, paths and custom functions
source("../../../source/main_packages.R")
source("../../../source/paths.R")
source("../../../source/functions.R")

# Load supplementary packages
packages <- c("RColorBrewer", "ggpubr", "cowplot", "Biostrings", "openxlsx", "kableExtra")
invisible(lapply(packages, require, character.only = TRUE))

Preparation

Tables preparation

Seqtab

# move to oligotyping directory
setwd(paste0(path_oligo,"/erwinia/oligotyping_Erwinia_sequences-c8-s1-a0.0-A0-M10"))

# load the matrix count table
matrix_count <- read.table("MATRIX-COUNT.txt", header = TRUE) %>% t()

# arrange it
colnames(matrix_count) <- matrix_count[1,]
matrix_count <- matrix_count[-1,]
matrix_count <- matrix_count %>% as.data.frame()

# print it
matrix_count %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
CTC1 CTC10 CTC11 CTC12 CTC13 CTC14 CTC15 CTC2 CTC3 CTC4 CTC5 CTC6 CTC7 CTC9 NP27 NP34 NP36 S146 S164 S165 S166 S167 S20 S21 S22 S24 S30 S31 S32 S33 S34 S35 S36 S37 S38 S44 S45 S46 S48 S49
AAGACTTA 2561 530 3566 332 1919 152 2825 2145 2268 147 11 3121 27 622 6 1 1 1 2 4 1 2 5 1 2 7 4069 4534 1857 3578 692 8 2513 377 19 26 46 7 32 8
TGAGTCGA 499 108 918 74 418 34 570 441 412 28 0 616 6 149 1 1 0 0 0 0 0 0 1 2 0 1 582 648 208 766 69 0 448 84 0 3 18 1 0 4
AAGACTTG 536 123 1017 94 402 45 588 442 386 27 3 760 7 123 0 0 0 0 0 0 0 0 0 0 0 0 15 13 2 3 2 0 4 1 0 0 0 1 0 0
TGAGTCGG 101 21 247 25 90 10 99 94 72 2 1 158 1 35 0 0 0 0 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 1 0 0 0 0 0 0
TGAGTCGT 28 4 63 9 36 3 35 24 19 1 0 41 2 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
TGAGCCGA 10 3 24 4 10 3 18 12 4 1 0 10 0 3 0 0 0 0 0 0 0 0 0 0 0 0 1 2 1 3 0 0 0 0 0 0 0 0 0 0
AAGACCTA 10 0 9 1 3 0 0 9 8 0 0 6 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 2 0 0 0 0 0 0 0
AGGACTTA 2 2 6 0 3 1 10 5 0 1 0 5 0 1 0 0 0 0 0 0 0 0 0 0 0 0 2 6 2 0 0 0 1 0 0 0 0 0 0 0
AAGGCTTA 7 1 4 0 5 0 5 2 3 0 0 6 0 5 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 2 0 0 3 1 0 0 0 0 0 0

Taxonomy

# move to oligotyping directory
setwd(paste0(path_oligo,"/erwinia/oligotyping_Erwinia_sequences-c8-s1-a0.0-A0-M10"))

# load the fasta table
fasta <- readDNAStringSet("OLIGO-REPRESENTATIVES.fasta")

# arrange it
fasta <- fasta %>% as.data.frame()
colnames(fasta) <- "seq"
fasta$oligotype <- rownames(fasta)
fasta <- fasta %>% dplyr::select(-c(seq))

# print it
fasta %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
oligotype
AAGACTTA AAGACTTA
TGAGTCGA TGAGTCGA
AAGACTTG AAGACTTG
TGAGTCGG TGAGTCGG
TGAGTCGT TGAGTCGT
TGAGCCGA TGAGCCGA
AAGACCTA AAGACCTA
AGGACTTA AGGACTTA
AAGGCTTA AAGGCTTA

Change oligotype name by oligotype / MED nodes in the matrix count

# Reference file 

## move to tsv directory
setwd(path_tsv)

## load the reference table
ref_oligo_med2 <- read.table("2B_REF_info_erwinia.tsv", sep="\t", header = TRUE)

## select only the 3 oligotypes of Erwinia
ref_oligo_med2 <- ref_oligo_med2[!is.na(ref_oligo_med2$oligotype),]

## change order of columns
ref_oligo_med2 <- ref_oligo_med2 %>% select(c(seq, oligotype, MED_node_frequency_size, OLIGO_oligotype_frequency_size))

## create a column with reference name (will be used in plots)
ref_oligo_med2$ref <- paste0("oligotype_", ref_oligo_med2$OLIGO_oligotype_frequency_size, " / node_", ref_oligo_med2$MED_node_frequency_size)

## create a copy of fasta 
fasta2 <- fasta

# Matrix count

## create an oligotype column in the matrix count
matrix_count$oligotype <- rownames(matrix_count)

## change order of columns
matrix_count <- matrix_count %>% dplyr::select(c(oligotype, everything()))

## merge the matrix count and the reference dataframe
matrix_count2 <- matrix_count %>% merge(ref_oligo_med2 %>% dplyr::select(-c(seq)), by="oligotype")

## change order of columns
matrix_count2 <- matrix_count2 %>% dplyr::select(c(oligotype, MED_node_frequency_size, OLIGO_oligotype_frequency_size, ref, everything()))

## change rownames
rownames(matrix_count2) <- matrix_count2$ref

## change order of columns
matrix_count2 <- matrix_count2 %>% dplyr::select(-c(oligotype, ref, MED_node_frequency_size, OLIGO_oligotype_frequency_size))

## print it
matrix_count2 %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
CTC1 CTC10 CTC11 CTC12 CTC13 CTC14 CTC15 CTC2 CTC3 CTC4 CTC5 CTC6 CTC7 CTC9 NP27 NP34 NP36 S146 S164 S165 S166 S167 S20 S21 S22 S24 S30 S31 S32 S33 S34 S35 S36 S37 S38 S44 S45 S46 S48 S49
oligotype_AAGACTTA (40) | size:38025 / node_N0798 (40) | size:38213 2561 530 3566 332 1919 152 2825 2145 2268 147 11 3121 27 622 6 1 1 1 2 4 1 2 5 1 2 7 4069 4534 1857 3578 692 8 2513 377 19 26 46 7 32 8
oligotype_AAGACTTG (22) | size:4594 / node_N0802 (22) | size:4627 536 123 1017 94 402 45 588 442 386 27 3 760 7 123 0 0 0 0 0 0 0 0 0 0 0 0 15 13 2 3 2 0 4 1 0 0 0 1 0 0
oligotype_TGAGTCGA (29) | size:7110 / node_N0311 (30) | size:8523 499 108 918 74 418 34 570 441 412 28 0 616 6 149 1 1 0 0 0 0 0 0 1 2 0 1 582 648 208 766 69 0 448 84 0 3 18 1 0 4 
## edit the fasta dataframe
fasta2 <- fasta2 %>% merge(ref_oligo_med2 %>% dplyr::select(-c(seq)), by="oligotype")
rownames(fasta2) <- fasta2$ref
fasta2 <- fasta2 %>% dplyr::select(-c(MED_node_frequency_size, OLIGO_oligotype_frequency_size, oligotype))

## print it
fasta2 %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
ref
oligotype_AAGACTTA (40) | size:38025 / node_N0798 (40) | size:38213 oligotype_AAGACTTA (40) | size:38025 / node_N0798 (40) | size:38213
oligotype_AAGACTTG (22) | size:4594 / node_N0802 (22) | size:4627 oligotype_AAGACTTG (22) | size:4594 / node_N0802 (22) | size:4627
oligotype_TGAGTCGA (29) | size:7110 / node_N0311 (30) | size:8523 oligotype_TGAGTCGA (29) | size:7110 / node_N0311 (30) | size:8523

Metadata

metadata <- read.csv(paste0(path_metadata,"/metadata_22_06_21.csv"), sep=";")
rownames(metadata) <- metadata$Sample

Phyloseq object with oligotypes

# convert matrix_count into matrix and numeric
matrix_count <- matrix_count2 %>% as.matrix()
class(matrix_count) <- "numeric"

# phyloseq elements
OTU = otu_table(as.matrix(matrix_count), taxa_are_rows =TRUE)
TAX = tax_table(as.matrix(fasta2))
SAM = sample_data(metadata)

# phyloseq object
ps <- phyloseq(OTU, TAX, SAM)
ps
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3 taxa and 40 samples ]
## sample_data() Sample Data:       [ 40 samples by 13 sample variables ]
## tax_table()   Taxonomy Table:    [ 3 taxa by 1 taxonomic ranks ]
compute_read_counts(ps)
## [1] 49729
# remove blanks
ps <- subset_samples(ps, Strain!="Blank")
ps <- check_ps(ps)
ps
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3 taxa and 40 samples ]
## sample_data() Sample Data:       [ 40 samples by 13 sample variables ]
## tax_table()   Taxonomy Table:    [ 3 taxa by 1 taxonomic ranks ]

Create new metadata with Percent

Load ps with all samples (for final plot)

setwd(path_rdata)
ps.filter <- readRDS("1D_MED_phyloseq_decontam.rds")
ps.filter <- check_ps(ps.filter)

Edit new metadata with Percent_erwinia

guide_italics <- guides(fill = guide_legend(label.theme = element_text(size = 16, face = "italic", colour = "Black", angle = 0)))

# add read depth in sample table of phyloseq object
sample_data(ps.filter)$Read_depth <- sample_sums(ps.filter)

# select Wolbachia
ps.erwinia <- ps.filter %>% subset_taxa(Genus=="Erwinia")

# add read depth of Wolbachia
sample_data(ps.filter)$Read_erwinia <- sample_sums(ps.erwinia)
sample_data(ps.filter) %>% colnames()
##  [1] "Sample"         "Well"           "Strain"         "Field"         
##  [5] "Country"        "Organ"          "Species"        "Individual"    
##  [9] "Individuals"    "Date"           "Run"            "Control"       
## [13] "Dna"            "Species_italic" "Strain_italic"  "Read_depth"    
## [17] "is.neg"         "Read_erwinia"
sample_data(ps.erwinia) %>% colnames()
##  [1] "Sample"         "Well"           "Strain"         "Field"         
##  [5] "Country"        "Organ"          "Species"        "Individual"    
##  [9] "Individuals"    "Date"           "Run"            "Control"       
## [13] "Dna"            "Species_italic" "Strain_italic"  "Read_depth"    
## [17] "is.neg"
# add percent of Wolbachia
sample_data(ps.filter)$Percent_erwinia <- sample_data(ps.filter)$Read_erwinia / sample_data(ps.filter)$Read_depth

# round the percent of Wolbachia at 2 decimals
sample_data(ps.filter)$Percent_erwinia <- sample_data(ps.filter)$Percent_erwinia %>% round(2)

# extract metadata table
test <- data.frame(sample_data(ps.filter))

# merge this metadata table with the other
new.metadata <- data.frame(sample_data(ps)) %>% merge(test %>% dplyr::select(c(Sample, Read_depth, Read_erwinia, Percent_erwinia)), by="Sample")
new.metadata <- test[new.metadata$Sample %in% sample_names(ps),]
rownames(new.metadata) <- new.metadata$Sample

# print it
new.metadata %>%
  kbl() %>%
  kable_paper("hover", full_width = F)
Sample Well Strain Field Country Organ Species Individual Individuals Date Run Control Dna Species_italic Strain_italic Read_depth is.neg Read_erwinia Percent_erwinia
CTC1 CTC1 G5 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 99 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 29779 FALSE 3769 0.13
CTC10 CTC10 D6 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 77 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 2609 FALSE 799 0.31
CTC11 CTC11 E6 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 99 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 13874 FALSE 5878 0.42
CTC12 CTC12 F6 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 87 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 1146 FALSE 542 0.47
CTC13 CTC13 G6 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 99 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 18035 FALSE 2898 0.16
CTC14 CTC14 H6 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 101 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 1708 FALSE 253 0.15
CTC15 CTC15 I6 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 99 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 23180 FALSE 4159 0.18
CTC2 CTC2 H5 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 99 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 30692 FALSE 3181 0.10
CTC3 CTC3 I5 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 99 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 39920 FALSE 3184 0.08
CTC4 CTC4 J5 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 29 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 2139 FALSE 207 0.10
CTC5 CTC5 K5 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 99 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 15789 FALSE 15 0.00
CTC6 CTC6 L5 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 99 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 19753 FALSE 4739 0.24
CTC9 CTC9 C6 Laboratory - Slab TC (Wolbachia -) Lab France Whole Culex quinquefasciatus N/A 0 09/03/2018 run2 True sample 67 italic(“Culex quinquefasciatus”) paste(“Laboratory - Slab TC (”, italic(“Wolbachia”), “-)”) 4980 FALSE 948 0.19
NP14 NP14 K4 Field - Guadeloupe Field Guadeloupe Ovary Aedes aegypti 1a 0 N/A run3 True sample 19 italic(“Aedes aegypti”) Field-Guadeloupe 7973 FALSE 0 0.00
NP2 NP2 K3 Field - Guadeloupe Field Guadeloupe Ovary Culex quinquefasciatus 1c 0 N/A run3 True sample 89 italic(“Culex quinquefasciatus”) Field-Guadeloupe 648335 FALSE 0 0.00
NP20 NP20 E5 Field - Guadeloupe Field Guadeloupe Ovary Aedes aegypti 3a 0 N/A run3 True sample 18 italic(“Aedes aegypti”) Field-Guadeloupe 136 FALSE 0 0.00
NP27 NP27 L5 Field - Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus 7c 0 N/A run3 True sample 30 italic(“Culex quinquefasciatus”) Field-Guadeloupe 1234 FALSE 7 0.01
NP29 NP29 B6 Field - Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus 9c 0 N/A run3 True sample 16 italic(“Culex quinquefasciatus”) Field-Guadeloupe 203 FALSE 0 0.00
NP30 NP30 C6 Field - Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus 10c 0 N/A run3 True sample 24 italic(“Culex quinquefasciatus”) Field-Guadeloupe 228 FALSE 0 0.00
NP34 NP34 G6 Field - Guadeloupe Field Guadeloupe Whole Culex quinquefasciatus 14c 0 N/A run3 True sample 20 italic(“Culex quinquefasciatus”) Field-Guadeloupe 95 FALSE 2 0.02
NP35 NP35 H6 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti 7a 0 N/A run3 True sample 84 italic(“Aedes aegypti”) Field-Guadeloupe 196532 FALSE 0 0.00
NP36 NP36 I6 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti 8a 0 N/A run3 True sample 27 italic(“Aedes aegypti”) Field-Guadeloupe 249 FALSE 1 0.00
NP37 NP37 J6 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti 9a 0 N/A run3 True sample 52 italic(“Aedes aegypti”) Field-Guadeloupe 419340 FALSE 0 0.00
NP38 NP38 K6 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti 10a 0 N/A run3 True sample 62 italic(“Aedes aegypti”) Field-Guadeloupe 282479 FALSE 0 0.00
NP39 NP39 L6 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti 11a 0 N/A run3 True sample 47 italic(“Aedes aegypti”) Field-Guadeloupe 218684 FALSE 0 0.00
NP41 NP41 B7 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti 13a 0 N/A run3 True sample 93 italic(“Aedes aegypti”) Field-Guadeloupe 247152 FALSE 0 0.00
NP42 NP42 C7 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti 14a 0 N/A run3 True sample 85 italic(“Aedes aegypti”) Field-Guadeloupe 185157 FALSE 0 0.00
NP43 NP43 D7 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti 15a 0 N/A run3 True sample 104 italic(“Aedes aegypti”) Field-Guadeloupe 239335 FALSE 0 0.00
NP44 NP44 E7 Field - Guadeloupe Field Guadeloupe Whole Aedes aegypti 16a 0 N/A run3 True sample 50 italic(“Aedes aegypti”) Field-Guadeloupe 156879 FALSE 0 0.00
NP5 NP5 B4 Field - Guadeloupe Field Guadeloupe Ovary Culex quinquefasciatus 2c 0 N/A run3 True sample 66 italic(“Culex quinquefasciatus”) Field-Guadeloupe 736159 FALSE 0 0.00
NP8 NP8 E4 Field - Guadeloupe Field Guadeloupe Ovary Culex quinquefasciatus 3c 0 N/A run3 True sample 92 italic(“Culex quinquefasciatus”) Field-Guadeloupe 334799 FALSE 0 0.00
S100 S100 K7 Field - Camping Europe Field France Ovary Culex pipiens GL1 1 30/05/2017 run1 True sample 105 italic(“Culex pipiens”) Field-Camping~Europe 52486 FALSE 0 0.00
S102 S102 A8 Field - Camping Europe Field France Ovary Culex pipiens GL2 2 30/05/2017 run1 True sample 14 italic(“Culex pipiens”) Field-Camping~Europe 3456 FALSE 0 0.00
S104 S104 C8 Field - Camping Europe Field France Ovary Culex pipiens GL5 5 30/05/2017 run1 True sample 54 italic(“Culex pipiens”) Field-Camping~Europe 52403 FALSE 0 0.00
S105 S105 D8 Field - Camping Europe Field France Ovary Culex pipiens GL6 6 30/05/2017 run1 True sample 103 italic(“Culex pipiens”) Field-Camping~Europe 55577 FALSE 0 0.00
S106 S106 E8 Field - Camping Europe Field France Ovary Culex pipiens GL7 7 30/05/2017 run1 True sample 96 italic(“Culex pipiens”) Field-Camping~Europe 33053 FALSE 0 0.00
S107 S107 F8 Field - Camping Europe Field France Ovary Culex pipiens GL8 8 30/05/2017 run1 True sample 65 italic(“Culex pipiens”) Field-Camping~Europe 52154 FALSE 0 0.00
S108 S108 G8 Field - Camping Europe Field France Ovary Culex pipiens GL9 9 30/05/2017 run1 True sample 64 italic(“Culex pipiens”) Field-Camping~Europe 55735 FALSE 0 0.00
S109 S109 H8 Field - Camping Europe Field France Ovary Culex pipiens GL10 10 30/05/2017 run1 True sample 55 italic(“Culex pipiens”) Field-Camping~Europe 59023 FALSE 0 0.00
S110 S110 I8 Field - Camping Europe Field France Ovary Culex pipiens GL11 0 30/05/2017 run1 True sample 57 italic(“Culex pipiens”) Field-Camping~Europe 57377 FALSE 0 0.00
S121 S121 H1 Field - Bosc Field France Ovary Culex pipiens J26 22 28/06/2017 run2 True sample 37 italic(“Culex pipiens”) Field-Bosc 20361 FALSE 0 0.00
S122 S122 I1 Field - Bosc Field France Ovary Culex pipiens J27 23 28/06/2017 run2 True sample 52 italic(“Culex pipiens”) Field-Bosc 9803 FALSE 0 0.00
S123 S123 J1 Field - Bosc Field France Ovary Culex pipiens J28 24 28/06/2017 run2 True sample 102 italic(“Culex pipiens”) Field-Bosc 20130 FALSE 0 0.00
S124 S124 K1 Field - Bosc Field France Ovary Culex pipiens J29 25 28/06/2017 run2 True sample 69 italic(“Culex pipiens”) Field-Bosc 18146 FALSE 0 0.00
S126 S126 K6 Field - Bosc Field France Ovary Culex pipiens J30 26 28/06/2017 run2 True sample 99 italic(“Culex pipiens”) Field-Bosc 15235 FALSE 0 0.00
S127 S127 B2 Field - Bosc Field France Ovary Culex pipiens J31 27 28/06/2017 run2 True sample 38 italic(“Culex pipiens”) Field-Bosc 24696 FALSE 0 0.00
S128 S128 C2 Field - Bosc Field France Ovary Culex pipiens J32 28 28/06/2017 run2 True sample 71 italic(“Culex pipiens”) Field-Bosc 16305 FALSE 0 0.00
S146 S146 I3 Laboratory - Lavar Lab France Ovary Culex pipiens MW52 29 29/08/2017 run2 True sample 72 italic(“Culex pipiens”) Laboratory-Lavar 25012 FALSE 1 0.00
S147 S147 J3 Laboratory - Lavar Lab France Ovary Culex pipiens MW53 30 29/08/2017 run2 True sample 56 italic(“Culex pipiens”) Laboratory-Lavar 25171 FALSE 0 0.00
S148 S148 K3 Laboratory - Lavar Lab France Ovary Culex pipiens MW54 31 29/08/2017 run2 True sample 91 italic(“Culex pipiens”) Laboratory-Lavar 14164 FALSE 0 0.00
S150 S150 A4 Laboratory - Lavar Lab France Ovary Culex pipiens MW55 32 29/08/2017 run2 True sample 43 italic(“Culex pipiens”) Laboratory-Lavar 15081 FALSE 0 0.00
S151 S151 B4 Laboratory - Lavar Lab France Ovary Culex pipiens MW56 33 29/08/2017 run2 True sample 78 italic(“Culex pipiens”) Laboratory-Lavar 22944 FALSE 0 0.00
S152 S152 C4 Laboratory - Lavar Lab France Ovary Culex pipiens MW57 34 29/08/2017 run2 True sample 79 italic(“Culex pipiens”) Laboratory-Lavar 15082 FALSE 0 0.00
S153 S153 D4 Laboratory - Lavar Lab France Ovary Culex pipiens MW58 35 29/08/2017 run2 True sample 97 italic(“Culex pipiens”) Laboratory-Lavar 17040 FALSE 0 0.00
S154 S154 E4 Laboratory - Lavar Lab France Ovary Culex pipiens MW59 36 29/08/2017 run2 True sample 76 italic(“Culex pipiens”) Laboratory-Lavar 9626 FALSE 0 0.00
S160 S160 K4 Laboratory - Lavar Lab France Ovary Culex pipiens MW60 37 29/08/2017 run2 True sample 99 italic(“Culex pipiens”) Laboratory-Lavar 72508 FALSE 0 0.00
S162 S162 B5 Laboratory - Lavar Lab France Ovary Culex pipiens MW61 38 29/08/2017 run2 True sample 90 italic(“Culex pipiens”) Laboratory-Lavar 25180 FALSE 0 0.00
S163 S163 L6 Laboratory - Lavar Lab France Ovary Culex pipiens MW62 39 30/08/2017 run2 True sample 96 italic(“Culex pipiens”) Laboratory-Lavar 12333 FALSE 0 0.00
S164 S164 C5 Laboratory - Lavar Lab France Ovary Culex pipiens MW63 40 30/08/2017 run2 True sample 74 italic(“Culex pipiens”) Laboratory-Lavar 22368 FALSE 2 0.00
S165 S165 D5 Laboratory - Lavar Lab France Ovary Culex pipiens MW64 41 30/08/2017 run2 True sample 98 italic(“Culex pipiens”) Laboratory-Lavar 17731 FALSE 4 0.00
S166 S166 E5 Field - Camping Europe Field France Ovary Culex pipiens GL4 4 30/05/2017 run2 True sample 53 italic(“Culex pipiens”) Field-Camping~Europe 13979 FALSE 1 0.00
S167 S167 F5 Field - Bosc Field France Ovary Culex pipiens J32 28 28/06/2017 run2 True sample 59 italic(“Culex pipiens”) Field-Bosc 14048 FALSE 2 0.00
S169 S169 B7 Field - Camping Europe Field France Ovary Culex pipiens 5 43 16/05/2017 run2 True sample 95 italic(“Culex pipiens”) Field-Camping~Europe 11553 FALSE 0 0.00
S170 S170 C7 Field - Camping Europe Field France Ovary Culex pipiens 6 44 16/05/2017 run2 True sample 94 italic(“Culex pipiens”) Field-Camping~Europe 8852 FALSE 0 0.00
S18 S18 A1 Laboratory - Lavar Lab France Whole Culex pipiens MW75 0 30/08/2017 run1 True sample 3 italic(“Culex pipiens”) Laboratory-Lavar 4290 FALSE 0 0.00
S19 S19 B1 Laboratory - Lavar Lab France Whole Culex pipiens MW65 0 30/08/2017 run1 True sample 51 italic(“Culex pipiens”) Laboratory-Lavar 44527 FALSE 0 0.00
S20 S20 C1 Laboratory - Lavar Lab France Whole Culex pipiens MW66 0 30/08/2017 run1 True sample 39 italic(“Culex pipiens”) Laboratory-Lavar 42864 FALSE 6 0.00
S21 S21 D1 Laboratory - Lavar Lab France Whole Culex pipiens MW67 0 30/08/2017 run1 True sample 36 italic(“Culex pipiens”) Laboratory-Lavar 33798 FALSE 3 0.00
S22 S22 E1 Laboratory - Lavar Lab France Whole Culex pipiens MW68 0 30/08/2017 run1 True sample 54 italic(“Culex pipiens”) Laboratory-Lavar 19044 FALSE 2 0.00
S23 S23 F1 Laboratory - Lavar Lab France Whole Culex pipiens MW69 0 30/08/2017 run1 True sample 48 italic(“Culex pipiens”) Laboratory-Lavar 38172 FALSE 0 0.00
S24 S24 G1 Laboratory - Lavar Lab France Whole Culex pipiens MW70 0 30/08/2017 run1 True sample 70 italic(“Culex pipiens”) Laboratory-Lavar 42355 FALSE 8 0.00
S25 S25 H1 Laboratory - Lavar Lab France Whole Culex pipiens MW71 0 30/08/2017 run1 True sample 51 italic(“Culex pipiens”) Laboratory-Lavar 47688 FALSE 0 0.00
S26 S26 I1 Laboratory - Lavar Lab France Whole Culex pipiens MW72 0 30/08/2017 run1 True sample 17 italic(“Culex pipiens”) Laboratory-Lavar 5394 FALSE 0 0.00
S27 S27 J1 Laboratory - Lavar Lab France Whole Culex pipiens MW73 0 30/08/2017 run1 True sample 34 italic(“Culex pipiens”) Laboratory-Lavar 24558 FALSE 0 0.00
S28 S28 A2 Laboratory - Lavar Lab France Whole Culex pipiens MW74 0 30/08/2017 run1 True sample 6 italic(“Culex pipiens”) Laboratory-Lavar 4503 FALSE 0 0.00
S30 S30 K1 Laboratory - Lavar Lab France Whole Culex pipiens MW1 0 23/08/2017 run1 True sample 82 italic(“Culex pipiens”) Laboratory-Lavar 25353 FALSE 4683 0.18
S31 S31 L1 Laboratory - Lavar Lab France Whole Culex pipiens MW2 0 23/08/2017 run1 True sample 44 italic(“Culex pipiens”) Laboratory-Lavar 20417 FALSE 5214 0.26
S32 S32 C2 Laboratory - Lavar Lab France Whole Culex pipiens MW3 0 23/08/2017 run1 True sample 22 italic(“Culex pipiens”) Laboratory-Lavar 12441 FALSE 2075 0.17
S33 S33 D2 Laboratory - Lavar Lab France Whole Culex pipiens MW4 0 23/08/2017 run1 True sample 26 italic(“Culex pipiens”) Laboratory-Lavar 33867 FALSE 4362 0.13
S34 S34 E2 Laboratory - Lavar Lab France Whole Culex pipiens MW5 0 23/08/2017 run1 True sample 31 italic(“Culex pipiens”) Laboratory-Lavar 9367 FALSE 764 0.08
S35 S35 F2 Laboratory - Lavar Lab France Whole Culex pipiens MW6 0 23/08/2017 run1 True sample 23 italic(“Culex pipiens”) Laboratory-Lavar 11663 FALSE 8 0.00
S36 S36 G2 Laboratory - Lavar Lab France Whole Culex pipiens MW7 0 23/08/2017 run1 True sample 79 italic(“Culex pipiens”) Laboratory-Lavar 33020 FALSE 2972 0.09
S37 S37 H2 Laboratory - Lavar Lab France Whole Culex pipiens MW8 0 23/08/2017 run1 True sample 88 italic(“Culex pipiens”) Laboratory-Lavar 18340 FALSE 464 0.03
S38 S38 I2 Laboratory - Lavar Lab France Whole Culex pipiens MW9 0 23/08/2017 run1 True sample 63 italic(“Culex pipiens”) Laboratory-Lavar 54790 FALSE 19 0.00
S39 S39 J2 Laboratory - Lavar Lab France Whole Culex pipiens MW10 0 23/08/2017 run1 True sample 75 italic(“Culex pipiens”) Laboratory-Lavar 36273 FALSE 0 0.00
S40 S40 K2 Laboratory - Lavar Lab France Whole Culex pipiens MW11 0 23/08/2017 run1 True sample 83 italic(“Culex pipiens”) Laboratory-Lavar 44448 FALSE 0 0.00
S42 S42 A3 Field - Camping Europe Field France Whole Culex pipiens GLE1 0 30/05/2017 run1 True sample 1 italic(“Culex pipiens”) Field-Camping~Europe 4107 FALSE 0 0.00
S43 S43 B3 Field - Camping Europe Field France Whole Culex pipiens GLE2 0 30/05/2017 run1 True sample 11 italic(“Culex pipiens”) Field-Camping~Europe 9279 FALSE 0 0.00
S44 S44 C3 Field - Camping Europe Field France Whole Culex pipiens GLE3 0 30/05/2017 run1 True sample 4 italic(“Culex pipiens”) Field-Camping~Europe 8026 FALSE 29 0.00
S45 S45 D3 Field - Camping Europe Field France Whole Culex pipiens GLE4 0 30/05/2017 run1 True sample 13 italic(“Culex pipiens”) Field-Camping~Europe 18150 FALSE 64 0.00
S47 S47 F3 Field - Camping Europe Field France Whole Culex pipiens GLE6 0 30/05/2017 run1 True sample 15 italic(“Culex pipiens”) Field-Camping~Europe 1951 FALSE 0 0.00
S48 S48 G3 Field - Camping Europe Field France Whole Culex pipiens GLE7 0 30/05/2017 run1 True sample 60 italic(“Culex pipiens”) Field-Camping~Europe 56738 FALSE 32 0.00
S49 S49 H3 Field - Bosc Field France Whole Culex pipiens E1 0 28/06/2017 run1 True sample 28 italic(“Culex pipiens”) Field-Bosc 33498 FALSE 12 0.00
S50 S50 I3 Field - Bosc Field France Whole Culex pipiens E2 0 28/06/2017 run1 True sample 25 italic(“Culex pipiens”) Field-Bosc 28481 FALSE 0 0.00
S51 S51 J3 Field - Bosc Field France Whole Culex pipiens E3 0 28/06/2017 run1 True sample 40 italic(“Culex pipiens”) Field-Bosc 61788 FALSE 0 0.00
S52 S52 K3 Field - Bosc Field France Whole Culex pipiens E4 0 28/06/2017 run1 True sample 21 italic(“Culex pipiens”) Field-Bosc 21553 FALSE 0 0.00
S55 S55 B4 Field - Bosc Field France Whole Culex pipiens E6 0 28/06/2017 run1 True sample 46 italic(“Culex pipiens”) Field-Bosc 50447 FALSE 0 0.00
S56 S56 C4 Field - Bosc Field France Whole Culex pipiens E7 0 28/06/2017 run1 True sample 61 italic(“Culex pipiens”) Field-Bosc 42609 FALSE 0 0.00
S57 S57 D4 Field - Bosc Field France Whole Culex pipiens E8 0 28/06/2017 run1 True sample 86 italic(“Culex pipiens”) Field-Bosc 49157 FALSE 0 0.00
S58 S58 E4 Field - Bosc Field France Whole Culex pipiens E9 0 28/06/2017 run1 True sample 32 italic(“Culex pipiens”) Field-Bosc 30357 FALSE 0 0.00
S59 S59 F4 Field - Bosc Field France Whole Culex pipiens E10 0 28/06/2017 run1 True sample 33 italic(“Culex pipiens”) Field-Bosc 32798 FALSE 0 0.00
S60 S60 G4 Field - Bosc Field France Whole Culex pipiens E11 0 28/06/2017 run1 True sample 45 italic(“Culex pipiens”) Field-Bosc 44485 FALSE 0 0.00
S61 S61 H4 Field - Bosc Field France Whole Culex pipiens E12 0 28/06/2017 run1 True sample 41 italic(“Culex pipiens”) Field-Bosc 49545 FALSE 0 0.00
S63 S63 J4 Field - Bosc Field France Whole Culex pipiens E14 0 28/06/2017 run1 True sample 100 italic(“Culex pipiens”) Field-Bosc 53444 FALSE 0 0.00
S64 S64 K4 Field - Bosc Field France Whole Culex pipiens E15 0 28/06/2017 run1 True sample 81 italic(“Culex pipiens”) Field-Bosc 47628 FALSE 0 0.00
S79 S79 B6 Field - Camping Europe Field France Ovary Culex pipiens J16 12 28/06/2017 run1 True sample 68 italic(“Culex pipiens”) Field-Camping~Europe 59755 FALSE 0 0.00
S80 S80 C6 Field - Camping Europe Field France Ovary Culex pipiens J17 13 28/06/2017 run1 True sample 83 italic(“Culex pipiens”) Field-Camping~Europe 52788 FALSE 0 0.00
S83 S83 F6 Field - Camping Europe Field France Ovary Culex pipiens J20 16 28/06/2017 run1 True sample 73 italic(“Culex pipiens”) Field-Camping~Europe 42272 FALSE 0 0.00
S84 S84 G6 Field - Camping Europe Field France Ovary Culex pipiens J21 17 28/06/2017 run1 True sample 49 italic(“Culex pipiens”) Field-Camping~Europe 56676 FALSE 0 0.00
S85 S85 H6 Field - Camping Europe Field France Ovary Culex pipiens J22 18 28/06/2017 run1 True sample 35 italic(“Culex pipiens”) Field-Camping~Europe 41690 FALSE 0 0.00
S86 S86 I6 Field - Camping Europe Field France Ovary Culex pipiens J23 19 28/06/2017 run1 True sample 80 italic(“Culex pipiens”) Field-Camping~Europe 61984 FALSE 0 0.00
S87 S87 J6 Field - Bosc Field France Ovary Culex pipiens J24 20 28/06/2017 run1 True sample 58 italic(“Culex pipiens”) Field-Bosc 65958 FALSE 0 0.00
S88 S88 K6 Field - Bosc Field France Ovary Culex pipiens J25 21 28/06/2017 run1 True sample 106 italic(“Culex pipiens”) Field-Bosc 53102 FALSE 0 0.00 
# replace metadata in the created phyloseq object
sample_data(ps) <- sample_data(new.metadata)

Taxonomic structure

Count

col <- brewer.pal(7, "Pastel2")

# reshape data for plot
test3 <- test %>% select(c(Sample, Species, Strain, Organ, Read_depth, Read_erwinia)) %>% reshape2::melt(id.vars=c("Sample", "Species", "Strain", "Organ"), vars=c("Read_depth", "Read_erwinia"))

count_whole <- test3[test3$Organ=="Whole",]
count_ovary <- test3[test3$Organ=="Ovary",]

make.italic <- function(x) as.expression(lapply(x, function(y) bquote(italic(.(y)))))

levels(count_whole$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(count_ovary$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(count_whole$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))

levels(count_ovary$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))


# plot
p_count1 <- ggplot(count_whole, aes(x = Sample, y = value, fill=variable))+ 
  geom_bar(position = "dodge", stat = "identity")+
  scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=12, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text=element_text(size=14), 
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
  facet_wrap(~Species+Strain+Organ, scales = "free_x", ncol=3, labeller=label_parsed)+
  labs(y="Sequence counts")+
  ylim(0, 900000)+
  geom_text(aes(label=value), position=position_dodge(width=1.1), width=0.25, size=4, hjust=-0.25, vjust=0.5, angle=90)+
  guides(fill=guide_legend(title="Read"))
## Warning: Ignoring unknown parameters: width
p_count2 <- ggplot(count_ovary, aes(x = Sample, y = value, fill=variable))+ 
  geom_bar(position = "dodge", stat = "identity")+
    scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=18, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text=element_text(size=14), 
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
 facet_wrap(~Species+Strain+Organ, scales = "free_x", ncol=3, labeller=label_parsed)+
  labs(y="Sequence counts")+
    ylim(0, 900000)+
  geom_text(aes(label=value), position=position_dodge(width=0.8), width=0.25, size=4, hjust=-0.25, vjust=0.5, angle=90)+
  guides(fill=guide_legend(title="Read"))
## Warning: Ignoring unknown parameters: width
# afficher plot
p_count1
## Warning: position_dodge requires non-overlapping x intervals
## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

p_count2

# panels
p_group <- plot_grid(p_count1+theme(legend.position="none"), 
          p_count2+theme(legend.position="none"), 
          nrow=2, 
          ncol=1)+
    draw_plot_label(c("B1", "B2"), c(0, 0), c(1, 0.5), size = 20)
## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals
legend_plot <- get_legend(p_count1 + theme(legend.position="bottom"))
## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals
p_counts <- plot_grid(p_group, legend_plot, nrow=2, ncol=1, rel_heights = c(1, .1))
p_counts

Whole (the most abundant nodes)

guide_italics <- guides(fill = guide_legend(label.theme = element_text(size = 16, face = "italic", colour = "Black", angle = 0),
                                            nrow=2, byrow=TRUE))

# select whole
ps.filter.whole <- subset_samples(ps, Organ=="Whole")
ps.filter.whole <- prune_taxa(taxa_sums(ps.filter.whole) >= 1, ps.filter.whole)
ps.filter.whole <- prune_samples(sample_sums(ps.filter.whole) >= 1, ps.filter.whole)
ps.filter.whole
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 3 taxa and 33 samples ]
## sample_data() Sample Data:       [ 33 samples by 19 sample variables ]
## tax_table()   Taxonomy Table:    [ 3 taxa by 1 taxonomic ranks ]
# data pour plot
#data_for_plot2 <- taxo_data_fast(ps.filter.whole, method = "abundance")
data_for_plot2 <- taxo_data(ps.filter.whole, top=15)
## Warning in psmelt(ps_global): The sample variables: 
## Sample
##  have been renamed to: 
## sample_Sample
## to avoid conflicts with special phyloseq plot attribute names.
paste0("\n15 MOST ABUNDANT GENUS: \n") %>% cat()
## 
## 15 MOST ABUNDANT GENUS:
paste0("\"", levels(data_for_plot2$Name), "\",\n") %>% cat()
## "oligotype_AAGACTTA (40) | size:38025 / node_N0798 (40) | size:38213.",
##  "oligotype_AAGACTTG (22) | size:4594 / node_N0802 (22) | size:4627.",
##  "oligotype_TGAGTCGA (29) | size:7110 / node_N0311 (30) | size:8523.",
data_for_plot2$Name <- data_for_plot2$Name %>% gsub(pattern = "node_", replacement ="" ) %>% as.factor()
data_for_plot2$Name <- as.factor(data_for_plot2$Name)

new_names <- c("oligotype_AAGACTTA (40) | size:38025 / N0798 (40) | size:38213.",
               "oligotype_TGAGTCGA (29) | size:7110 / N0311 (30) | size:8523.",
               "oligotype_AAGACTTG (22) | size:4594 / N0802 (22) | size:4627.",
               "Other.")

data_for_plot2$Name <- factor(data_for_plot2$Name, levels = new_names)

col_add <- brewer.pal(8, "Accent")

col <- c("oligotype_AAGACTTA (40) | size:38025 / N0798 (40) | size:38213."="#FFFFCF",
         "oligotype_TGAGTCGA (29) | size:7110 / N0311 (30) | size:8523."="#FFE352",
         "oligotype_AAGACTTG (22) | size:4594 / N0802 (22) | size:4627."="#F5F61B",
         "Other."="#A0A0A0")

levels(data_for_plot2$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
               "Culex pipiens"=make.italic("Culex pipiens"),
               "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))

levels(data_for_plot2$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))

#data_for_plot2 <- data_for_plot2 %>% na.omit()

p2 <- ggplot(data_for_plot2, aes(x = Sample, y = Relative_Abundance, fill = Name, species=Species, organ=Organ, Strain=Strain))+ 
  geom_bar(position = "stack", stat = "identity")+
  scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=18, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text = element_text(size=14),
        #legend.key.height = unit(1, 'cm'),
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
  facet_wrap(~Species+Strain+Organ, scales = "free", ncol=3, labeller=label_parsed)+
  labs(x="Sample", y="Relative abundance", fill="Oligotype / MED node")

p2

Ovary (the most abundant nodes)

# select ovary
ps.filter.ovary <- subset_samples(ps, Organ=="Ovary")
ps.filter.ovary <- prune_taxa(taxa_sums(ps.filter.ovary) >= 1, ps.filter.ovary)
ps.filter.ovary <- prune_samples(sample_sums(ps.filter.ovary) >= 1, ps.filter.ovary)
ps.filter.ovary
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1 taxa and 5 samples ]
## sample_data() Sample Data:       [ 5 samples by 19 sample variables ]
## tax_table()   Taxonomy Table:    [ 1 taxa by 1 taxonomic ranks ]
# data pour plot
data_for_plot3 <- taxo_data(ps.filter.ovary, top=15)
## Warning in psmelt(ps_global): The sample variables: 
## Sample
##  have been renamed to: 
## sample_Sample
## to avoid conflicts with special phyloseq plot attribute names.
paste0("\n15 MOST ABUNDANT GENUS: \n") %>% cat()
## 
## 15 MOST ABUNDANT GENUS:
paste0("\"", levels(data_for_plot3$Name), "\",\n") %>% cat()
## "oligotype_AAGACTTA (40) | size:38025 / node_N0798 (40) | size:38213.",
data_for_plot3$Name <- data_for_plot3$Name %>% gsub(pattern = "node_", replacement ="" ) %>% as.factor()
data_for_plot3$Name <- as.factor(data_for_plot3$Name)

levels(data_for_plot3$Species)= c("Culex pipiens"=make.italic("Culex pipiens"))

levels(data_for_plot3$Strain) <- c("Bosc", "Camping~Europe",  "Lavar~(lab)")

p3 <- ggplot(data_for_plot3, aes(x = Sample, y = Relative_Abundance, fill = Name, species=Species, organ=Organ, Strain=Strain))+ 
  geom_bar(position = "stack", stat = "identity")+
  scale_fill_manual(values = col)+
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, size=18, hjust=1, vjust=0.5)) +
  ggtitle("") + 
  guide_italics+
  theme(legend.title = element_text(size = 20), 
        legend.position="bottom",
        legend.text = element_text(size=14),
        #legend.key.height = unit(1, 'cm'),
        panel.spacing.y=unit(1, "lines"), 
        panel.spacing.x=unit(0.8, "lines"),
        panel.spacing=unit(0,"lines"),
        strip.background=element_rect(color="grey30", fill="grey90"),
        strip.text.x = element_text(size = 16),
        panel.border=element_rect(color="grey90"),
        axis.ticks.x=element_blank(),
        axis.text.y = element_text(size=18)) +
  facet_wrap(~Species+Strain+Organ, scales = "free", ncol=3, labeller=label_parsed)+
  labs(x="Sample", y="Relative abundance", fill="Oligotype / MED node")

p3

# ps.filter.ovary@otu_table # only 1 oligotype
# 
# a <- data.frame(ps.filter.ovary@sam_data)
# b <- data.frame(ps.filter.ovary@tax_table)
# count_table <- data.frame(ps.filter.ovary@otu_table)
# 
# count_table %>%
#   kbl() %>%
#   kable_paper("hover", full_width = F)
# 
# a$Relative_Abundance <- 1.0000000000
# a$Name <- paste0(b$ref,".")
# 
# levels(a$Species)= c("Aedes aegypti"=make.italic("Aedes aegypti"),
#                "Culex pipiens"=make.italic("Culex pipiens"),
#                "Culex quinquefasciatus"=make.italic("Culex quinquefasciatus"))
# 
# levels(a$Strain) <- c("Bosc", "Camping~Europe", "Guadeloupe", "Lavar~(lab)", expression(paste(italic("Wolbachia"), "- (Slab TC)")))
# 
# p3 <- ggplot(a, aes(x = Sample, y = Relative_Abundance, fill = Name, species=Species, organ=Organ, Strain=Strain))+
#   geom_bar(position = "stack", stat = "identity")+
#   scale_fill_manual(values = "#FFFFCF")+
#   theme_bw() +
#   theme(axis.text.x = element_text(angle = 90)) +
#   ggtitle("") +
#   guide_italics+
#   theme(legend.title = element_text(size = 20), 
#         legend.position="bottom",
#         legend.text = element_text(size=14),
#         #legend.key.height = unit(1, 'cm'),
#         panel.spacing.y=unit(1, "lines"), 
#         panel.spacing.x=unit(0.8, "lines"),
#         panel.spacing=unit(0,"lines"),
#         strip.background=element_rect(color="grey30", fill="grey90"),
#         strip.text.x = element_text(size = 16),
#         panel.border=element_rect(color="grey90"),
#         axis.ticks.x=element_blank(),
#         axis.text.y = element_text(size=18)) +
#  facet_wrap(~Species+Strain+Organ, scales = "free", ncol=3, labeller=label_parsed)+
#   labs(x="Sample", y="Relative abundance", fill="Oligotype / MED Node")
# 
# p3

Panels taxonomy of whole / ovary

legend_plot <- get_legend(p2 + theme(legend.position="bottom"))

# panels
p_group <- plot_grid(p2+theme(legend.position="none"), 
          p3+theme(legend.position="none"), 
          nrow=2, 
          ncol=1)+
    draw_plot_label(c("A1", "A2"), c(0, 0), c(1, 0.5), size = 20)

p_taxo <- plot_grid(p_group, legend_plot, nrow=2, rel_heights = c(1, .1))
p_taxo

Save taxonomic plot

setwd(path_plot)

tiff("2De_OLIGO_counts_erwinia.tiff", units="in", width=20, height=18, res=300)
p_counts
dev.off()
## quartz_off_screen 
##                 2
tiff("2De_OLIGO_taxonomic_erwinia_whole.tiff", units="in", width=16, height=12, res=300)
p2
dev.off()
## quartz_off_screen 
##                 2
tiff("2De_OLIGO_taxonomic_erwinia_ovary.tiff", units="in", width=18, height=14, res=300)
p3
dev.off()
## quartz_off_screen 
##                 2
tiff("2De_OLIGO_taxonomic_erwinia.tiff", units="in", width=18, height=16, res=300)
p_taxo
dev.off()
## quartz_off_screen 
##                 2
png("2De_OLIGO_counts_erwinia_big.png", units="in", width=20, height=18, res=300)
p_counts
dev.off()
## quartz_off_screen 
##                 2
png("2De_OLIGO_counts_erwinia_small.png", units="in", width=18, height=14, res=300)
p_counts
dev.off()
## quartz_off_screen 
##                 2
png("2De_OLIGO_taxonomic_erwinia_whole.png", units="in", width=16, height=12, res=300)
p2
dev.off()
## quartz_off_screen 
##                 2
png("2De_OLIGO_taxonomic_erwinia_ovary.png", units="in", width=18, height=14, res=300)
p3
dev.off()
## quartz_off_screen 
##                 2
png("2De_OLIGO_taxonomic_erwinia_big.png", units="in", width=18, height=18, res=300)
p_taxo
dev.off()
## quartz_off_screen 
##                 2
png("2De_OLIGO_taxonomic_erwinia_small.png", units="in", width=18, height=14, res=300)
p_taxo
dev.off()
## quartz_off_screen 
##                 2

Make main plot

setwd(paste0(path_oligo,"/erwinia/oligotyping_Erwinia_sequences-c8-s1-a0.0-A0-M10/HTML-OUTPUT"))

img <- magick::image_read("entropy.png")
p_entropy <- magick::image_ggplot(img, interpolate = TRUE)
p_entropy+ theme(plot.margin = unit(c(-7,-2.5,-7,-0.5), "cm"))

p_entropy+ theme(plot.margin=unit(c(-7,-2,-12,-5), "mm"))

aligned <- plot_grid(p_taxo, 
                     p_counts, 
                     align="hv")

aligned

p_entropy2 <- plot_grid(p_entropy, nrow=1)+
  draw_plot_label(c("C"), c(0), c(1), size=20, hjust=-0.5)

p_entropy2

t_plot <- plot_grid(aligned, 
                    p_entropy2,
                    nrow=2, 
                    ncol=1, 
                    scale=1,
                    rel_heights=c(2,1))

t_plot

setwd(path_plot)

tiff("2De_OLIGO_main_erwinia.tiff", width=36, height=36, res=300, units="in")
t_plot
dev.off()
## quartz_off_screen 
##                 2
png("2De_OLIGO_main_erwinia.png", width=36, height=36, res=300, units="in")
t_plot
dev.off()
## quartz_off_screen 
##                 2